%%% Matteo Ciccarelli and Benoit Mojon (c).  Replicate results of Table 3
%%% in the paper

clear all;

% load a panel of GDP for G7 countries
data = xlsread('cpi_DATA');
% ordering of the series
%series = {CPIQU2	CPIQUS	CPIQCA	CPIQUK	CPIQJP	CPIQDE	CPIQFR	CPIQIT	CPIQAT	CPIQBE	CPIQDK	CPIQFI	CPIQGR	CPIQIE	CPIQLU	CPIQPT	CPIQES	CPIQSE	CPIQNL	CPIQAU	CPIQNZ	CPIQNO	CPIQCH	CPIOECD


nco = 24;
Y = data(:,1:nco); 
n = size(Y,2);
dly = [log(Y(5:end,:))-log(Y(1:end-4,:))]*100;

filter = bpf(dly,6,32,12);

x = filter(1+12:end-12,2:n-1);
%x = filter(9:178,2:n-1);
[T,N] = size(x);

%%% principal components
OPTS.disp=0;
Mx = mean(x);
Wx = (std(x));

%%% center data
x_std = (x-kron(ones(T,1),Mx))*diag(1./Wx);

[v,d] = eigs(cov(x_std),1,'lm',OPTS);
F_pc = x_std*v*d^-.5;
chi_pc = x_std*v*v'*diag(Wx) + kron(ones(size(x,1),1),Mx);

%%% factor analysis
r = 1; %% number of factors

%[chi_zz,F_zz,C_zz,R_zz] = fa_zz(x,r); %% static em etimation
[chi_em,F_em,C_em,R_em] = FA_em(x,r); %% static zig-zag estimation
%[chi_dyn,F_dyn,A_dyn, C_dyn, Q_dyn, R_dyn] = dyn_em(x,r,0,r,1); %% dynamic factor analysis q = r = 1; p =1; s = 0
%[chi_dyn4,F_dyn4,A_dyn4, C_dyn4, Q_dyn4, R_dyn4] = dyn_em(x,r,1,r,1); %% dynamic factor analysis q = r = 1; p =2; s = 3

% cross-sectional averages
avg_x = mean(x_std,2); avg_x = avg_x./std(avg_x);
chi_mean = avg_x*inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx) + kron(ones(size(x,1),1),Mx);
C_ave = (inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx))';
%%% compute sign to rotate estimated factors
%temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
%temp = corrcoef(avg_x,F_zz(:,1)); s_zz=sign(temp(2,1));
%temp = corrcoef(avg_x,F_em(:,1)); s_em=sign(temp(2,1));
%temp = corrcoef(avg_x,F_dyn(:,1)); s_dyn=sign(temp(2,1));

% temp = corrcoef([avg_x,F_dyn]); s_dyn = sign(temp(1,2:end));  
%temp = corrcoef([avg_x,F_zz]); s_zz = sign(temp(1,2:end));  
temp = corrcoef([avg_x,F_em]); s_em = sign(temp(1,2:end));  
% temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
if s_em(1,1)<0
    F_em(:,1) = s_em(1,1)*F_em(:,1);
end

% if s_zz(1,1)<0
%     F_zz(:,1) = s_zz(1,1)*F_zz(:,1);
% end
    
    vari_em1=(C_em(:,1).^2)*cov(F_em(:,1))./diag(cov(x));
%    vari_zz=(C_zz(:,1).^2)*cov(F_zz(:,1))./diag(cov(x));
%     vari_em2=(C_em(:,2).^2)*cov(F_em(:,2))./diag(cov(x));
%     disp('Variance decomposition')
%     disp([ vari_em1 vari_em2 ])
 %   disp([ vari_em1 ])
 
% if r==1
%     F_dyn = s_dyn*F_dyn ;
%     F_em = s_em*F_em ;
%     vari_ave=(C_ave(:,1).^2)*cov(avg_x(:,1))./diag(cov(x));
%     vari_em=(C_em(:,1).^2)*cov(F_em(:,1))./diag(cov(x));
%     vari_dyn=(C_dyn(:,1).^2)*cov(F_dyn(:,1))./diag(cov(x));
%     disp([vari_ave vari_em vari_dyn])
% end 
%     
% if r==2
%     vari_dyn1=(C_dyn(:,1).^2)*cov(F_dyn(:,1))./diag(cov(x));
%     vari_dyn2=(C_dyn(:,2).^2)*cov(F_dyn(:,2))./diag(cov(x));
%     disp('Variance decomposition')
%     disp([ vari_dyn1 vari_dyn2 ])
% end
%
x = filter(1+12:end-12,1:n-1);
%x = filter(9:178,1:n-1);
[T,N] = size(x);
oecd_fact = (filter(1+12:end-12,nco)-kron(ones(T,1),mean(filter(1+12:end-12,nco))))*diag(1./std(filter(1+12:end-12,nco)));
%oecd_fact = (filter(1+12:178,nco)-kron(ones(T,1),mean(filter(9:178,nco))))*diag(1./std(filter(9:178,nco)));

regr = [ones(T,1) oecd_fact];

for i= 1:N
    par=inv(regr'*regr)*regr'*x(:,i);           %%% ols coeff
    u=x(:,i)-regr*par;                    %%% residuals
    m=size(regr); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_oecd(i,1)=1-(u'*u)/(center(x(:,i))'*center(x(:,i)));              %%% R2
    r2_oecd(i,2)=1-(1-r2_oecd(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

regr_av = [ones(T,1) avg_x];

for i= 1:N
    par=inv(regr_av'*regr_av)*regr_av'*x(:,i);           %%% ols coeff
    u=x(:,i)-regr_av*par;                    %%% residuals
    m=size(regr_av); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_av(i,1)=1-(u'*u)/(center(x(:,i))'*center(x(:,i)));              %%% R2
    r2_av(i,2)=1-(1-r2_av(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

regr_fa = [ones(T,1) F_em(:,1)];
par=inv(regr_fa'*regr_fa)*regr_fa'*x(:,1);           %%% ols coeff
u=x(:,1)-regr_fa*par;                    %%% residuals
m=size(regr_fa); Ti=m(1); k=m(2);
esu=u'*u/(Ti-k);
r2_eu1(1,1)=1-(u'*u)/(center(x(:,1))'*center(x(:,1)));              %%% R2
r2_eu1(1,2)=1-(1-r2_eu1(1,1))*(Ti-1)/(Ti-k);        %%% adjusted R2

% Table 3
sh_var = [r2_av(:,1) r2_oecd(:,1) ];
sh_fa(2:23,1) = vari_em1(1:22,1);
sh_fa(1,1)=r2_eu1(1,1);
disp('table 3')
table3_nosort = [sh_var sh_fa]

